Systems are known in the art for computing numerical solutions of linear systems of equations. Also known is the preconditioning of the coefficient matrix of such a system.
Computer aided design (CAD) systems have become accepted universally as important engineering design tools for a variety of applications in mechanical, civil, and electrical engineering. Information regarding the shape and physical properties of a device or structure are captured and stored in the system in the form of a geometric model. The behavior of the physical system is described in mathematical terms, via the application of known laws of physics, as a problem whose solution reveals a property of the system that it is desired to evaluate. In practice, a region is divided into elements, often tens of thousands or millions of elements, over which the properties are considered to vary in a known manner and to which classical principles can be applied to yield a system of linear algebraic equations.
The computational cost associated with evaluating the physical system is dominated by the solution of this system of equations. It is therefore of paramount technical importance to improve the computational efficiency with which the numerical solution to these equations can be generated. For this reason, iterative techniques are typically the method of choice for solving such linear systems due to their speed advantage over direct inversion methods. These methods, however, are strongly dependent of the condition number, i.e., the size of the span of the spectrum of the resulting matrix. It is known that large condition numbers reduce the convergence of the system dramatically.
In order to overcome the convergence rate difficulties, matrix preconditioners are used to reduce the condition number. Solving linear systems using conjugate gradient (CG) iterative methods, for example, with preconditioning has been well documented and is summarized in “Numerical Recipes”, W. Press et al., Cambridge Press (1992). The most common preconditioning method is known as incomplete LU, and is described in detail in “Iterative Solution Methods”, O. Axelsson, Cambridge Press (1994). To summarize, in a triangular (LU) factorization of a matrix A, the matrix is decomposed into a lower triangular matrix L and an upper triangular matrix U, thereby expressing the matrix as a product LU=A. If A is a sparse matrix, meaning the matrix has a low ratio of nonzero elements, the ratio of nonzero elements in the resulting matrices L and U is increased. The generation of nonzero elements is ordinarily called “fill-in”. In the incomplete LU, the fill-in is ignored by approximating it to zero. Alternatively, fill-in can be selectively ignored based on its relative magnitude. Modifications to the incomplete LU in generating preconditioners for sparse coefficient matrices are described in U.S. Pat. No. 5,136,538, “Preconditioned Conjugate Gradient System”, N. K. Karmardar et al., and in commonly assigned U.S. Pat. No. 5,754,181, “Computer Aided Design System”, V. Amdursky et al.
The application of incomplete LU is limited to instances where the condition number is small, which is especially the case in the solution of dense matrices, such as those produced when attempting to solve integral equations. In order to apply an incomplete LU on a dense preconditioner, the preconditioner must first be “sparsed”. For example, setting to zero all matrix elements whose magnitude is below a certain threshold, or setting to zero all elements in the preconditioner representing interactions between elements in the physical system separated by distances larger then a certain threshold. All known choices of sparsing the preconditioner, including the two mentioned above, when combined with incomplete, or even complete, LU factorization produce poor results difficulty in general linear systems with large conditioner numbers, a general and robust preconditioning method is needed. Prior to this invention, this need was not adequately met.
A particular problem arises during the electrical analysis of devices that are small with respect to the wavelength, such as electronic packaging structures, chips, interconnects, and printed circuit boards. The analysis of such structures results in linear systems of equations to be solved, which are desired to be solved iteratively by a fast, low storage requirement solver such as the Fast Multipole Method or the pre-corrected FFT method. The use of an iterative solver requires that the condition number of the linear system (the ratio of the largest to smallest eigenvalue) be as close to one as possible. When the condition number is too large a preconditioner is used to condition the system and reduce the condition number. However, for such structures as those mentioned above, the use of common rooftop basis functions (also known as Rao-Wilton Glisson or RWG functions) produces linear systems with large condition numbers such that no known preconditoner has been available to properly condition the system.
As was noted above, the most common attempt at a preconditioner is one based on geometrically close interactions. To form such a preconditioner, one sets to zero all interactions in the linear system that are geometrically separated by a defined distance. Such a preconditioner will only be effective, however, for structures that are large compared to the wavelength when the common rooftop or RWG basis functions are used. One existing solution to this problem uses less desirable loop-tree or loop-star basis functions which result in a better conditioned linear system. A disadvantage of this approach, however, is that the fast solvers mentioned above have difficulties dealing with the often large loops produced by the change in the basis functions. This is due to the inaccurate calculation of near interactions by the fast solver, for which large loops have many. Reference can be made to “Fast and efficient algorithms in electromagneties”, by Weng Chew et al., Artech House 2001, for descriptions of fast solvers and loop-tree or loop-star basis functions.